#include "base.h"

int main(int argc,char **argv)
{

//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@///

   //////////////////
   // Input   ///////
   //////////////////
   char dir[1024];
   char gammaprint[8];

   int  set_p[4];
   char setpprint[16];
   char Oa2print[16];
   int iprop;
   int beta100;

   char momentum_order[1024];
 ////// Parsing command line

   for(int i=1; i<argc; i++) {
            if(strncmp(argv[i],"-Nclust",7)==0)  sscanf(argv[i],"-Nclust=%d",&Nclust);
            if(strncmp(argv[i],"-dir",4)==0)  sscanf(argv[i],"-dir=%s",dir);
            if(strncmp(argv[i],"-T",2)==0)  sscanf(argv[i],"-T=%d",&T);
            if(strncmp(argv[i],"-L",2)==0)  sscanf(argv[i],"-L=%d",&L);
	    if(strncmp(argv[i],"-gamma_to_prom",14)==0)  sscanf(argv[i],"-gamma_to_prom=%s",gammaprint);    

	    if(strncmp(argv[i],"-beta100",8)==0)  sscanf(argv[i],"-beta100=%d",&beta100);
            if(strncmp(argv[i],"-set_p",6)==0)  sscanf(argv[i],"-set_p=%d,%d,%d,%d",&set_p[3], &set_p[2], &set_p[1], &set_p[0]);

            if(strncmp(argv[i],"-iprop",6)==0)  sscanf(argv[i],"-iprop=%d",&iprop);
            if(strncmp(argv[i],"-Oa2",4)==0)  sscanf(argv[i],"-Oa2=%s",Oa2print);

            if(strncmp(argv[i],"-momentum_order",15)==0)  sscanf(argv[i],"-momentum_order=%s",momentum_order);

  }

   /////////////////////////
  //// Momenta set
  //////////////////////////////

  whichProp(iprop,setpprint);
  char momentafile[1024]; sprintf(momentafile,"../INPUT/mom_%d.dat",beta100);
  read_momenta_data(momentafile);
  double *nn=(double*)malloc(sizeof(double)*(Nmom));
  compute_nn(momentum_order,nn);
  Momentum_Magro(momentum_order,set_p);
  // equivalent momenta
  Nmom_eq=Number_of_Equivalent_Momenta(SinP2);
  Average_over_Equivalent_Momenta(SinP2,SinP2,SinP2_eq);

  //////////////////////////////////
  // IO files
  ///////////////////////////////

  char file[1024];
  char dir_qav[1024];
  char pm[8];

//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//



   double ***Z=malloc(sizeof(double**)*(3));
   for (int ipm=0; ipm<3; ipm++) Z[ipm]=malloc(sizeof(double*)*(Nmom_eq));
   for (int ipm=0; ipm<3; ipm++) for (int imom=0; imom<Nmom_eq; imom++) Z[ipm][imom]=malloc(sizeof(double)*(Nclust+1));

   double *aux_clust=malloc(sizeof(double)*(Nclust+1));

    // read 
    for (int ipm=0; ipm<2; ipm++)
    {
    whichpm(ipm,pm);
    sprintf(file,"Z%s_%s",gammaprint,Oa2print);
    sprintf(dir_qav,"%s%s/%s/QAV/",dir,pm,setpprint);
    for (int imom=0; imom<Nmom_eq; imom++)
    {
		read_input_onemomenta(dir_qav,file,aux_clust,imom);
    		for (int iclust=0; iclust<Nclust+1; iclust++) Z[ipm][imom][iclust]=aux_clust[iclust];
    }
    }

   // prom
   for (int imom=0; imom<Nmom_eq; imom++) for (int iclust=0; iclust<Nclust+1; iclust++) Z[2][imom][iclust]=(Z[0][imom][iclust]+Z[1][imom][iclust])/2.;
    


   // print output
   sprintf(dir_qav,"%sa/%s/QAV/",dir,setpprint);
   printf("Writting output in %s\n",dir_qav);
         sprintf(file,"Z%s_%s",gammaprint,Oa2print);
         for (int imom=0; imom<Nmom_eq; imom++) {
                for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust[iclust]= Z[2][imom][iclust];
                 write_output_onemomenta_real(dir_qav,file,aux_clust,SinP2_eq,imom);
          }
          sprintf(file,"%s/Z_ALLP/Z%s_%s.dat",dir_qav,gammaprint,Oa2print);
          FILE *fout=fopen(file,"w");
          if(fout == (FILE *)NULL) {
              fprintf(stderr,"Unable to open file %s for writing\n",file);
              return EXIT_FAILURE;
          }
         for (int imom=0; imom<Nmom_eq; imom++) {
                for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust[iclust]=Z[2][imom][iclust];
                fprintf(fout,"%8.10f %8.10f %8.10f  \n",SinP2_eq[imom],Z[2][imom][Nclust],error_jack_real(aux_clust,Nclust));
        }
      fclose(fout);


   return 0;

}
